/*************************************************************
** TRUBIT1: Truncated bivariate probit model
** Monte Carlo Comparison of Trubit with Boolean Probit
** (c) 2003 by Sanford C. Gordon
**
**
** Department of Politics
** new York University
** Purpose: Determine sample properties of partially observed probit with Alastair's trubit model.
** output file: trubit1.out
*************************************************************/

rndseed 33333;

closeall; 
cls; 

/* Call to the Maxlik library */
library maxlik;
/* Resets global variables to default values */
maxset;

/* Global Parameter Settings */

_max_algorithm = 2;
	                 	@ Optimization Method           @
                       	@ 1 = Steepest Descent          @
                       	@ 2 = BFGS                      @
                       	@ 3 = DFP                       @
                       	@ 4 = Newton-Raphson            @
                       	@ 5 = BHHH                      @
                       	@ 6 = PRCG                      @

_max_CovPar = 1;
						@ Covariance Matrix Method		@
						@ 0 = Not Computed			    @
						@ 1 = Numerical Hessian			@
						@ 2 = Quasi-maximum likelihood  @

_max_LineSearch = 2;
						@ Line Search Method			@
						@ 1 = unit step length			@
						@ 2 = cubic or quadratic		@
						@ 3 = step halving			    @
						@ 4 = Brent's method			@
						@ 5 = BHHH method				@
/*__output=0;*/
/*Create symbolic infinity*/
infty=2^10000;

proc throbiti(parms,dta);
    local obs,y,y1,y2,x1,x2,beta1,beta2,tau1,tau2,idx1,idx2,p0,py1,py2,pya,llike;
    /* Data */
    obs = rows(dta);
    y   = dta[.,1];
    y1  = dta[.,2];
    y2  = dta[.,3];
    x1  = ones(obs,1)~dta[.,4:5];
    x2  = ones(obs,1)~dta[.,4]~dta[.,6];
    /* parameters */
    beta1 = parms[1:3];
    beta2 = parms[4:6];
    tau1  = exp(parms[7]);
    tau2  = exp(parms[8]);
    idx1= x1*beta1;
    idx2=x2*beta2;
    p0 = cdfn(-idx1).*cdfn(-idx2);
    py1 = cdfn(idx1-tau1).*cdfn(tau2-idx2);
    py2 = cdfn(idx2-tau2).*cdfn(tau1-idx1);
    pya = 1-p0-py1-py2;
    llike = (1-y).*ln(p0)+y.*y1.*ln(py1)+y.*y2.*ln(py2)+y.*(1-y1).*(1-y2).*ln(pya);
    retp(llike);
endp;

proc boolprob(parms,dta);
    local y,llike,obs,x1,x2,beta1,beta2,p1,p2,yesprob,noprob;
    obs=rows(dta);
    y=dta[.,1];
    x1  = ones(obs,1)~dta[.,4:5];
    x2  = ones(obs,1)~dta[.,4]~dta[.,6];
    beta1=parms[1:3];
    beta2=parms[4:6];
    p1 = cdfn(x1*beta1);
    p2 = cdfn(x2*beta2);
    yesprob = 1-(1-p1).*(1-p2);
    noprob  = 1-yesprob;
    llike = y.*ln(yesprob)+(1-y).*ln(noprob);
    retp(llike);
endp;

proc trubiti(parms,dta);
    local y,y1,y2,llike,obs,x1,x2,beta1,beta2,p1,p2,yesproba,yesprob1,yesprob2,noprob;
    obs=rows(dta);
    y=dta[.,1];
    y1=dta[.,2];
    y2=dta[.,3];
    x1 = ones(obs,1)~dta[.,4:5];
    x2 = ones(obs,1)~dta[.,4]~dta[.,6];
    beta1 = parms[1:3];
    beta2 = parms[4:6];
    p1 = cdfn(x1*beta1);
    p2 = cdfn(x2*beta2);
    noprob = (1-p1).*(1-p2);
    yesprob1 = p1;
    yesprob2 = p2;
    yesproba = 1-noprob;
    llike = (1-y).*ln(noprob) + y.*(1-y1).*(1-y2).*ln(yesproba) + y.*y1.*(1-y2).*ln(yesprob1)
            +y.*(1-y1).*y2.*ln(yesprob2);
    retp(llike);
endp;

/* Parameter Values Fixed for all simulations */

beta1=1|1|1;
beta2=1|1|1;
tau1 = 2;
tau2 = 2;
rho=0;
mciter=1000;

/**************************************************************************************
** First Monte Carlo Routine: Compare trubit and Boolean probit, for sample of 100  **
**************************************************************************************/
/* Generate Xs, which are fixed in repeated samples */
N=100;
xcommon = rndn(N,1);
x12 = rndn(N,1)-1.7;
x22 = rndn(N,1)-1.7;
x1 = ones(N,1)~xcommon~x12;
x2 = ones(N,1)~xcommon~x22;

i=1; do until i >mciter;
cls;
print "N=" n ", i=" i;
/* Generate observed outcomes */
    epsilon1 = rndn(N,1);
    epsilon2 = rndn(N,1);
    zstar1  = x1*beta1+epsilon1;
    zstar2  = x2*beta2+epsilon2;
    y = maxc((zstar1~zstar2)').>0;
    reveal = rndu(n,1).>0.75;
/* Note how the trubit data generating process differs from throbit */
   
    y1 = (zstar1.>0).*reveal;
    y2 = (zstar2.>0).*reveal;
    both = y1.*y2;
    y1 = y1.*(1-both);
    y2 = y2.*(1-both);
 
    let label = y y1 y2 xcommon x12 x22;
    create f0= "c:/gauss50/research/throbit/trubit1" with ^label,0,8;
    call writer(f0,y~y1~y2~xcommon~x12~x22);
    call close(f0);

    dta = "c:/gauss50/research/throbit/trubit1"; /* Tells Maxlik where to look for the data */
    vars = getname(dta);
    
    /* First, run Boolean probit */
    _max_parnames = "beta10"|"beta1c"|"beta12"|"beta20"|"beta2c"|"beta22";
    theta0bp=0|0|0|0|0|0;
    {thetabp,fbp,gbp,covbp,retbp}=maxprt(maxlik(dta,vars,&boolprob,theta0bp));

    /* Now, run trubit */
    _max_parnames = "beta10"|"beta1c"|"beta12"|"beta20"|"beta2c"|"beta22";
    theta0tr=0|0|0|0|0|0;
    {thetatr,ftr,gtr,covtr,rettr}=maxprt(maxlik(dta,vars,&trubiti,theta0tr));
    
    /* Save the results */
    miss6 = {. . . . . .};
    
    if i == 1;
        meany = meanc(y~y1~y2)';
        /* Results for Boolean probit */
        resultbp = thetabp';
        convrgbp = retbp;
        llikebp = fbp*N;
        truthbp = ones(6,1);
        msematbp = (thetabp-truthbp)*(thetabp-truthbp)';
        if ismiss(covbp) == 0;
            serrbp = sqrt(diag(covbp))'; 
        else;
            serrbp = miss6;
        endif;
        /* Results for trubit */
        resulttr = thetatr';
        convrgtr = rettr;
        lliketr = ftr*N;
        truthtr = truthbp;
        msemattr = (thetatr-truthtr)*(thetatr-truthtr)';
        if ismiss(covtr) == 0;
            serrtr = sqrt(diag(covtr))'; 
        else;
            serrtr = miss6;
        endif;
    else;
        meany = meany|meanc(y~y1~y2)';
        resultbp = resultbp|thetabp';
        convrgbp = convrgbp|retbp;
        llikebp = llikebp|fbp*n;
        truthbp = ones(6,1);
        msematbp = msematbp+(thetabp-truthbp)*(thetabp-truthbp)';
        if ismiss(covbp)==0;
            serrbp=serrbp|sqrt(diag(covbp))';
        else;
            serrbp = serrbp|miss6;
        endif;

        resulttr = resulttr|thetatr';
        convrgtr = convrgtr|rettr;
        lliketr = lliketr|ftr*N;
        truthtr = truthbp;
        msemattr = msemattr+(thetatr-truthtr)*(thetatr-truthtr)';
        if ismiss(covtr)==0;
            serrtr=serrtr|sqrt(diag(covtr))';
        else;
            serrtr = serrtr|miss6;
        endif;
    endif;
    rmbp = meanc(resultbp);
    rmtr = meanc(resulttr);
    i=i+1;
endo;    

sebp100 = meanc(packr(serrbp));
setr100 = meanc(packr(serrtr));

rmsebp = sqrt(diag(msematbp))/mciter;
rmsetr = sqrt(diag(msemattr))/mciter;

/* Print table */ 
amper4= "&"$|"&"$|"&"$|"&";
slashes4  = "\\\\"$|"\\\\"$|"\\\\"$|"\\\\";
caption = "Truth"|"Mean Boolean"|"Mean Trubit"|"Rel. Eff.";
meanbp = meanc(resultbp);
meantr = meanc(resulttr);
releff = rmsebp./(rmsetr[1:6]);
b10 = truthbp[1]|meanbp[1]|meantr[1]|releff[1];
b1c = truthbp[2]|meanbp[2]|meantr[2]|releff[2];
b12 = truthbp[3]|meanbp[3]|meantr[3]|releff[3];
b20 = truthbp[4]|meanbp[4]|meantr[4]|releff[4]; 
b2c = truthbp[5]|meanbp[5]|meantr[5]|releff[5];
b22 = truthbp[6]|meanbp[6]|meantr[6]|releff[6];

b10s = "" $+ftocv(b10,0,2);
b1cs = "" $+ftocv(b1c,0,2);
b12s = "" $+ftocv(b12,0,2);
b20s = "" $+ftocv(b20,0,2);
b2cs = "" $+ftocv(b2c,0,2);
b22s = "" $+ftocv(b22,0,2);

resmat = caption$~amper4$~b10s$~amper4$~b1cs$~amper4$~b12s$~amper4$~b20s$~amper4$~b2cs$~amper4$~b22s$~slashes4;
format 200,2;
output file = "c:/gauss50/research/throbit/mc1tr.out" reset;
print"\\begin{table}";
print"\\caption{\\label{mc1} Monte Carlo Comparison of Boolean Probit and Trubit}";
print"\\begin{center}";
print"\\begin{tabular}{lcccccc}";
print"\\hline \\hline";
print"&$\\beta_{10}$&$\\beta_{1c}$&$\\beta_{13}$&$\\beta_{20}$&$\\beta_{2c}$&$\\beta_{23}$ \\\\ \\hline";

print "\\multicolumn{1}{l}{} &  \\multicolumn{6}{c}{\\emph{N=100}} \\\\ \\hline";
resmat;
print "\\hline";
print "\\multicolumn{7}{l}{Proportion successful convergence, Boolean Probit: " ""$+ftocv(rows(packr(serrbp))/rows(serrbp),0,2)$+"}\\\\";
print "\\multicolumn{7}{l}{Proportion successful convergence, Trubit: " ""$+ftocv(rows(packr(serrtr))/rows(serrtr),0,2)$+"}\\\\";
output off;

output file = "c:/gauss50/research/throbit/diagtr.out" reset;
print "I have finished with N=100, moving on to N=500";
output off;



/**************************************************************************************
** Second Monte Carlo Routine: Compare throbit and Boolean probit, for sample of 500 **
**************************************************************************************/
/* Generate Xs, which are fixed in repeated samples */
N=500;
xcommon = rndn(N,1);
x12 = rndn(N,1)-1.7;
x22 = rndn(N,1)-1.7;
x1 = ones(N,1)~xcommon~x12;
x2 = ones(N,1)~xcommon~x22;

i=1; do until i >mciter;
cls;
print "N=" n ", i=" i;
/* Generate observed outcomes */
    epsilon1 = rndn(N,1);
    epsilon2 = rndn(N,1);
    zstar1  = x1*beta1+epsilon1;
    zstar2  = x2*beta2+epsilon2;
    y = maxc((zstar1~zstar2)').>0;
    reveal = rndu(n,1).>0.75;
/* Note how the trubit data generating process differs from throbit */
    y1 = (zstar1.>0).*reveal;
    y2 = (zstar2.>0).*reveal;
    both = y1.*y2;
    y1 = y1.*(1-both);
    y2 = y2.*(1-both);
    let label = y y1 y2 xcommon x12 x22;
    create f0= "c:/gauss50/research/throbit/trubit1" with ^label,0,8;
    call writer(f0,y~y1~y2~xcommon~x12~x22);
    call close(f0);

    dta = "c:/gauss50/research/throbit/trubit1"; /* Tells Maxlik where to look for the data */
    vars = getname(dta);
    
    /* First, run Boolean probit */
    _max_parnames = "beta10"|"beta1c"|"beta12"|"beta20"|"beta2c"|"beta22";
    theta0bp=0|0|0|0|0|0;
    {thetabp,fbp,gbp,covbp,retbp}=maxprt(maxlik(dta,vars,&boolprob,theta0bp));

    /* Now, run trubit */
    _max_parnames = "beta10"|"beta1c"|"beta12"|"beta20"|"beta2c"|"beta22";
    theta0tr=0|0|0|0|0|0;
    {thetatr,ftr,gtr,covtr,rettr}=maxprt(maxlik(dta,vars,&trubiti,theta0tr));
    
    /* Save the results */
    miss6 = {. . . . . .};
    
    if i == 1;
        meany = meanc(y~y1~y2)';
        /* Results for Boolean probit */
        resultbp = thetabp';
        convrgbp = retbp;
        llikebp = fbp*N;
        truthbp = ones(6,1);
        msematbp = (thetabp-truthbp)*(thetabp-truthbp)';
        if ismiss(covbp) == 0;
            serrbp = sqrt(diag(covbp))'; 
        else;
            serrbp = miss6;
        endif;
        /* Results for trubit */
        resulttr = thetatr';
        convrgtr = rettr;
        lliketr = ftr*N;
        truthtr = truthbp;
        msemattr = (thetatr-truthtr)*(thetatr-truthtr)';
        if ismiss(covtr) == 0;
            serrtr = sqrt(diag(covtr))'; 
        else;
            serrtr = miss6;
        endif;
    else;
        meany = meany|meanc(y~y1~y2)';
        resultbp = resultbp|thetabp';
        convrgbp = convrgbp|retbp;
        llikebp = llikebp|fbp*n;
        truthbp = ones(6,1);
        msematbp = msematbp+(thetabp-truthbp)*(thetabp-truthbp)';
        if ismiss(covbp)==0;
            serrbp=serrbp|sqrt(diag(covbp))';
        else;
            serrbp = serrbp|miss6;
        endif;

        resulttr = resulttr|thetatr';
        convrgtr = convrgtr|rettr;
        lliketr = lliketr|ftr*N;
        truthtr = truthbp;
        msemattr = msemattr+(thetatr-truthtr)*(thetatr-truthtr)';
        if ismiss(covtr)==0;
            serrtr=serrtr|sqrt(diag(covtr))';
        else;
            serrtr = serrtr|miss6;
        endif;
    endif;
    i=i+1;
endo;    

sebp500 = meanc(packr(serrbp));
setr500 = meanc(packr(serrtr));

rmsebp = sqrt(diag(msematbp))/mciter;
rmsetr = sqrt(diag(msemattr))/mciter;

/* Print table continuation */ 
amper4= "&"$|"&"$|"&"$|"&";
slashes4  = "\\\\"$|"\\\\"$|"\\\\"$|"\\\\";
caption = "Truth"|"Mean Boolean"|"Mean Trubit"|"Rel. Eff.";
meanbp = meanc(resultbp);
meantr = meanc(resulttr);
releff = rmsebp./(rmsetr[1:6]);
b10 = truthbp[1]|meanbp[1]|meantr[1]|releff[1];
b1c = truthbp[2]|meanbp[2]|meantr[2]|releff[2];
b12 = truthbp[3]|meanbp[3]|meantr[3]|releff[3];
b20 = truthbp[4]|meanbp[4]|meantr[4]|releff[4]; 
b2c = truthbp[5]|meanbp[5]|meantr[5]|releff[5];
b22 = truthbp[6]|meanbp[6]|meantr[6]|releff[6];

b10s = "" $+ftocv(b10,0,2);
b1cs = "" $+ftocv(b1c,0,2);
b12s = "" $+ftocv(b12,0,2);
b20s = "" $+ftocv(b20,0,2);
b2cs = "" $+ftocv(b2c,0,2);
b22s = "" $+ftocv(b22,0,2);

resmat = caption$~amper4$~b10s$~amper4$~b1cs$~amper4$~b12s$~amper4$~b20s$~amper4$~b2cs$~amper4$~b22s$~slashes4;
output file = "c:/gauss50/research/throbit/mc1tr.out" on;
print "\\hline";
print "\\multicolumn{1}{l}{} &  \\multicolumn{6}{c}{\\emph{N=500}} \\\\ \\hline";
resmat;
print "\\hline";
print "\\multicolumn{7}{l}{Proportion successful convergence, Boolean Probit: " ""$+ftocv(rows(packr(serrbp))/rows(serrbp),0,2)$+"}\\\\";
print "\\multicolumn{7}{l}{Proportion successful convergence, Trubit: " ""$+ftocv(rows(packr(serrtr))/rows(serrtr),0,2)$+"}\\\\";
output off;

output file = "c:/gauss50/research/throbit/diagtr.out" on;
print "I have finished with N=500, moving on to N=1000";
output off;

/**************************************************************************************
** Third Monte Carlo Routine: Compare throbit and Boolean probit, for sample of 1000 **
**************************************************************************************/
/* Generate Xs, which are fixed in repeated samples */
N=1000;
xcommon = rndn(N,1);
x12 = rndn(N,1)-1.7;
x22 = rndn(N,1)-1.7;
x1 = ones(N,1)~xcommon~x12;
x2 = ones(N,1)~xcommon~x22;

i=1; do until i >mciter;
cls;
print "N=" n ", i=" i;
/* Generate observed outcomes */
    epsilon1 = rndn(N,1);
    epsilon2 = rndn(N,1);
    zstar1  = x1*beta1+epsilon1;
    zstar2  = x2*beta2+epsilon2;
    y = maxc((zstar1~zstar2)').>0;
    reveal = rndu(n,1).>0.75;
/* Note how the trubit data generating process differs from throbit */
    y1 = (zstar1.>0).*reveal;
    y2 = (zstar2.>0).*reveal;
    both = y1.*y2;
    y1 = y1.*(1-both);
    y2 = y2.*(1-both);
    let label = y y1 y2 xcommon x12 x22;
    create f0= "c:/gauss50/research/throbit/trubit1" with ^label,0,8;
    call writer(f0,y~y1~y2~xcommon~x12~x22);
    call close(f0);

    dta = "c:/gauss50/research/throbit/trubit1"; /* Tells Maxlik where to look for the data */
    vars = getname(dta);
    
    /* First, run Boolean probit */
    _max_parnames = "beta10"|"beta1c"|"beta12"|"beta20"|"beta2c"|"beta22";
    theta0bp=0|0|0|0|0|0;
    {thetabp,fbp,gbp,covbp,retbp}=maxprt(maxlik(dta,vars,&boolprob,theta0bp));

    /* Now, run trubit */
    _max_parnames = "beta10"|"beta1c"|"beta12"|"beta20"|"beta2c"|"beta22";
    theta0tr=0|0|0|0|0|0;
    {thetatr,ftr,gtr,covtr,rettr}=maxprt(maxlik(dta,vars,&trubiti,theta0tr));
    
    /* Save the results */
    miss6 = {. . . . . .};
    
    if i == 1;
        meany = meanc(y~y1~y2)';
        /* Results for Boolean probit */
        resultbp = thetabp';
        convrgbp = retbp;
        llikebp = fbp*N;
        truthbp = ones(6,1);
        msematbp = (thetabp-truthbp)*(thetabp-truthbp)';
        if ismiss(covbp) == 0;
            serrbp = sqrt(diag(covbp))'; 
        else;
            serrbp = miss6;
        endif;
        /* Results for trubit */
        resulttr = thetatr';
        convrgtr = rettr;
        lliketr = ftr*N;
        truthtr = truthbp;
        msemattr = (thetatr-truthtr)*(thetatr-truthtr)';
        if ismiss(covtr) == 0;
            serrtr = sqrt(diag(covtr))'; 
        else;
            serrtr = miss6;
        endif;
    else;
        meany = meany|meanc(y~y1~y2)';
        resultbp = resultbp|thetabp';
        convrgbp = convrgbp|retbp;
        llikebp = llikebp|fbp*n;
        truthbp = ones(6,1);
        msematbp = msematbp+(thetabp-truthbp)*(thetabp-truthbp)';
        if ismiss(covbp)==0;
            serrbp=serrbp|sqrt(diag(covbp))';
        else;
            serrbp = serrbp|miss6;
        endif;

        resulttr = resulttr|thetatr';
        convrgtr = convrgtr|rettr;
        lliketr = lliketr|ftr*N;
        truthtr = truthbp;
        msemattr = msemattr+(thetatr-truthtr)*(thetatr-truthtr)';
        if ismiss(covtr)==0;
            serrtr=serrtr|sqrt(diag(covtr))';
        else;
            serrtr = serrtr|miss6;
        endif;
    endif;
    i=i+1;
endo;    

sebp1000 = meanc(packr(serrbp));
setr1000 = meanc(packr(serrth));

rmsebp = sqrt(diag(msematbp))/mciter;
rmsetr = sqrt(diag(msemattr))/mciter;

/* Print table continuation */ 
amper4= "&"$|"&"$|"&"$|"&";
slashes4  = "\\\\"$|"\\\\"$|"\\\\"$|"\\\\";
caption = "Truth"|"Mean Boolean"|"Mean Trubit"|"Rel. Eff.";
meanbp = meanc(resultbp);
meantr = meanc(resulttr);
releff = rmsebp./(rmsetr[1:6]);
b10 = truthbp[1]|meanbp[1]|meantr[1]|releff[1];
b1c = truthbp[2]|meanbp[2]|meantr[2]|releff[2];
b12 = truthbp[3]|meanbp[3]|meantr[3]|releff[3];
b20 = truthbp[4]|meanbp[4]|meantr[4]|releff[4]; 
b2c = truthbp[5]|meanbp[5]|meantr[5]|releff[5];
b22 = truthbp[6]|meanbp[6]|meantr[6]|releff[6];

b10s = "" $+ftocv(b10,0,2);
b1cs = "" $+ftocv(b1c,0,2);
b12s = "" $+ftocv(b12,0,2);
b20s = "" $+ftocv(b20,0,2);
b2cs = "" $+ftocv(b2c,0,2);
b22s = "" $+ftocv(b22,0,2);

resmat = caption$~amper4$~b10s$~amper4$~b1cs$~amper4$~b12s$~amper4$~b20s$~amper4$~b2cs$~amper4$~b22s$~slashes4;
output file = "c:/gauss50/research/throbit/mc1tr.out" on;
print "\\hline";
print "\\multicolumn{1}{l}{} &  \\multicolumn{6}{c}{\\emph{N=1,000}} \\\\ \\hline";
resmat;
print "\\hline";
print "\\multicolumn{7}{l}{Proportion successful convergence, Boolean Probit: " ""$+ftocv(rows(packr(serrbp))/rows(serrbp),0,2)$+"}\\\\";
print "\\multicolumn{7}{l}{Proportion successful convergence, Trubit: " ""$+ftocv(rows(packr(serrtr))/rows(serrtr),0,2)$+"}\\\\";
output off;

output file = "c:/gauss50/research/throbit/diagtr.out" on;
print "I have finished with N=1000, moving on to N=5000";
output off;

/**************************************************************************************
** Fourth Monte Carlo Routine: Compare throbit and Boolean probit, for sample of 5000 **
**************************************************************************************/
/* Generate Xs, which are fixed in repeated samples */
N=5000;
xcommon = rndn(N,1);
x12 = rndn(N,1)-1.7;
x22 = rndn(N,1)-1.7;
x1 = ones(N,1)~xcommon~x12;
x2 = ones(N,1)~xcommon~x22;

i=1; do until i >mciter;
cls;
print "N=" n ", i=" i;
/* Generate observed outcomes */
    epsilon1 = rndn(N,1);
    epsilon2 = rndn(N,1);
    zstar1  = x1*beta1+epsilon1;
    zstar2  = x2*beta2+epsilon2;
    y = maxc((zstar1~zstar2)').>0;
    reveal = rndu(n,1).>0.75;
/* Note how the trubit data generating process differs from throbit */
    y1 = (zstar1.>0).*reveal;
    y2 = (zstar2.>0).*reveal;
    both = y1.*y2;
    y1 = y1.*(1-both);
    y2 = y2.*(1-both);
    let label = y y1 y2 xcommon x12 x22;
    create f0= "c:/gauss50/research/throbit/trubit1" with ^label,0,8;
    call writer(f0,y~y1~y2~xcommon~x12~x22);
    call close(f0);

    dta = "c:/gauss50/research/throbit/trubit1"; /* Tells Maxlik where to look for the data */
    vars = getname(dta);
    
    /* First, run Boolean probit */
    _max_parnames = "beta10"|"beta1c"|"beta12"|"beta20"|"beta2c"|"beta22";
    theta0bp=0|0|0|0|0|0;
    {thetabp,fbp,gbp,covbp,retbp}=maxprt(maxlik(dta,vars,&boolprob,theta0bp));

    /* Now, run trubit */
    _max_parnames = "beta10"|"beta1c"|"beta12"|"beta20"|"beta2c"|"beta22";
    theta0tr=0|0|0|0|0|0;
    {thetatr,ftr,gtr,covtr,rettr}=maxprt(maxlik(dta,vars,&trubiti,theta0tr));
    
    /* Save the results */
    miss6 = {. . . . . .};
    
    if i == 1;
        meany = meanc(y~y1~y2)';
        /* Results for Boolean probit */
        resultbp = thetabp';
        convrgbp = retbp;
        llikebp = fbp*N;
        truthbp = ones(6,1);
        msematbp = (thetabp-truthbp)*(thetabp-truthbp)';
        if ismiss(covbp) == 0;
            serrbp = sqrt(diag(covbp))'; 
        else;
            serrbp = miss6;
        endif;
        /* Results for trubit */
        resulttr = thetatr';
        convrgtr = rettr;
        lliketr = ftr*N;
        truthtr = truthbp;
        msemattr = (thetatr-truthtr)*(thetatr-truthtr)';
        if ismiss(covtr) == 0;
            serrtr = sqrt(diag(covtr))'; 
        else;
            serrtr = miss6;
        endif;
    else;
        meany = meany|meanc(y~y1~y2)';
        resultbp = resultbp|thetabp';
        convrgbp = convrgbp|retbp;
        llikebp = llikebp|fbp*n;
        truthbp = ones(6,1);
        msematbp = msematbp+(thetabp-truthbp)*(thetabp-truthbp)';
        if ismiss(covbp)==0;
            serrbp=serrbp|sqrt(diag(covbp))';
        else;
            serrbp = serrbp|miss6;
        endif;

        resulttr = resulttr|thetatr';
        convrgtr = convrgtr|rettr;
        lliketr = lliketr|ftr*N;
        truthtr = truthbp;
        msemattr = msemattr+(thetatr-truthtr)*(thetatr-truthtr)';
        if ismiss(covtr)==0;
            serrtr=serrtr|sqrt(diag(covtr))';
        else;
            serrtr = serrtr|miss6;
        endif;
    endif;
    i=i+1;
endo;    

sebp5000 = meanc(packr(serrbp));
setr5000 = meanc(packr(serrtr));

rmsebp = sqrt(diag(msematbp))/mciter;
rmsetr = sqrt(diag(msemattr))/mciter;

/* Print table continuation */ 
amper4= "&"$|"&"$|"&"$|"&";
slashes4  = "\\\\"$|"\\\\"$|"\\\\"$|"\\\\";
caption = "Truth"|"Mean Boolean"|"Mean Trubit"|"Rel. Eff.";
meanbp = meanc(resultbp);
meantr = meanc(resulttr);
releff = rmsebp./(rmsetr[1:6]);
b10 = truthbp[1]|meanbp[1]|meantr[1]|releff[1];
b1c = truthbp[2]|meanbp[2]|meantr[2]|releff[2];
b12 = truthbp[3]|meanbp[3]|meantr[3]|releff[3];
b20 = truthbp[4]|meanbp[4]|meantr[4]|releff[4]; 
b2c = truthbp[5]|meanbp[5]|meantr[5]|releff[5];
b22 = truthbp[6]|meanbp[6]|meantr[6]|releff[6];

b10s = "" $+ftocv(b10,0,2);
b1cs = "" $+ftocv(b1c,0,2);
b12s = "" $+ftocv(b12,0,2);
b20s = "" $+ftocv(b20,0,2);
b2cs = "" $+ftocv(b2c,0,2);
b22s = "" $+ftocv(b22,0,2);

N = 100000;
xcommon = rndn(N,1);
x12 = rndn(N,1)-1.7;
x22 = rndn(N,1)-1.7;
x1 = ones(N,1)~xcommon~x12;
x2 = ones(N,1)~xcommon~x22;

/* Generate observed outcomes */
    epsilon1 = rndn(N,1);
    epsilon2 = rndn(N,1);
    zstar1  = x1*beta1+epsilon1;
    zstar2  = x2*beta2+epsilon2;
    y = maxc((zstar1~zstar2)').>0;
    reveal = rndu(n,1).>0.75;
/* Note how the trubit data generating process differs from throbit */
    y1 = (zstar1.>0).*reveal;
    y2 = (zstar2.>0).*reveal;
    both = y1.*y2;
    y1 = y1.*(1-both);
    y2 = y2.*(1-both);

yprop  = meanc(y-y1-y2);
y1prop = meanc(y1);
y2prop = meanc(y2);

resmat = caption$~amper4$~b10s$~amper4$~b1cs$~amper4$~b12s$~amper4$~b20s$~amper4$~b2cs$~amper4$~b22s$~slashes4;
output file = "c:/gauss50/research/throbit/mc1tr.out" on;
print "\\hline";
print "\\multicolumn{1}{l}{} &  \\multicolumn{6}{c}{\\emph{N=5,000}} \\\\ \\hline";
resmat;
print "\\hline";
print "\\multicolumn{7}{l}{Proportion successful convergence, Boolean Probit: " ""$+ftocv(rows(packr(serrbp))/rows(serrbp),0,2)$+"}\\\\";
print "\\multicolumn{7}{l}{Proportion successful convergence, Trubit: " ""$+ftocv(rows(packr(serrtr))/rows(serrtr),0,2)$+"}\\\\";
print " \\hline \\hline ";
print "\\multicolumn{7}{l}{Approximate proportion ambiguous success: " ""$+ftocv(yprop,0,2)$+"}\\\\";
print "\\multicolumn{7}{l}{Approximate proportion success for first reason: " ""$+ftocv(y1prop,0,2)$+"}\\\\";
print "\\multicolumn{7}{l}{Approximate proportion success for second reason: " ""$+ftocv(y2prop,0,2)$+"}\\\\";
print "\\end{tabular}";
print "\\end{center}";
print "\\end{table}";
output off;

output file = "c:/gauss50/research/throbit/diagtr.out" on;
print "I have finished with N=5000.";
print "Mean standard errors, N=100";
print "Boolean probit:" sebp100';
print "trubit:" setr100';

print "Mean standard errors, N=500";
print "Boolean probit:" sebp500';
print "trubit:" setr500';

print "Mean standard errors, N=1000";
print "Boolean probit:" sebp1000';
print "trubit:" setr1000';

print "Mean standard errors, N=5000";
print "Boolean probit:" sebp5000';
print "trubit:" setr5000';
output off;





